/**********************************************************************************

	ran1, gasdev, sortd, indexxd


	Recipes from Numerical Recipes in C

 (C) Copr. 1986-92 Numerical Recipes Software


**********************************************************************************/


#include <math.h>
#include "nrutil.h"
#include "c_recipes.h"



/**********************************************************************************
   Returns a uniform random deviate between 0.0 and 1.0
**********************************************************************************/
#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

float ran1(long *idum)
/*
“Minimal” random number generator of Park and Miller with Bays-Durham shuffle and added
safeguards. Returns a uniform random deviate between 0.0 and 1.0 (exclusive of the endpoint
values). Call with idum a negative integer to initialize; thereafter, do not alter idum between
successive deviates in a sequence. RNMX should approximate the largest floating value that is
less than 1.
*/
{
	int j;
	long k;
	static long iy=0;
	static long iv[NTAB];
	float temp;

	if (*idum <= 0 || !iy) {			    /* Initialize */
		if (-(*idum) < 1) *idum=1;		  /* Be sure to prevent idum = 0 */
		else *idum = -(*idum);
		for (j=NTAB+7;j>=0;j--) {		    /* Load the shuffle table (after 8 warm-ups */
			k=(*idum)/IQ;
			*idum=IA*(*idum-k*IQ)-IR*k;
			if (*idum < 0) *idum += IM;
			if (j < NTAB) iv[j] = *idum;
		}
		iy=iv[0];
	}
	k=(*idum)/IQ;					            /* Start here when not initializing. */
	*idum=IA*(*idum-k*IQ)-IR*k;			        /* Compute idum=(IA*idum) % IM without over- */
	if (*idum < 0) *idum += IM;			        /*  flows by Schrage's method. */
	j=iy/NDIV;					                /* Will be in the range 0..NTAB-1. */
	iy=iv[j];					                /* Output previously stored value and refill the */
	iv[j] = *idum;					            /*  shuffle table. */
	if ((temp=AM*iy) > RNMX) return RNMX;		/* Because users don't expect endpoint values. */
	else return temp;
}

#undef IA
#undef IM
#undef AM
#undef IQ
#undef IR
#undef NTAB
#undef NDIV
#undef EPS
#undef RNMX





/**********************************************************************************
 Sorts an array arr[1..n] into ascending numerical order
**********************************************************************************/
#define NRANSI
#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
#define M 7
#define NSTACK 50
/* Here M is the size of subarrays sorted by straight insertion and NSTACK is the required */
/* auxiliary storage.*/


void sortd(unsigned long n, double arr[])
/* Sorts an array arr[1..n] into ascending numerical order using the Quicksort algorithm. n is */
/* input; arr is replaced on output by its sorted rearrangement. */

{
	unsigned long i,ir=n,j,k,l=1;
	int jstack=0,*istack;
	double a,temp;

	istack=ivector(1,NSTACK);
	for (;;) {						/* Insertion sort when subarray small enough. */
		if (ir-l < M) {
			for (j=l+1;j<=ir;j++) {
				a=arr[j];
				for (i=j-1;i>=1;i--) {
					if (arr[i] <= a) break;
					arr[i+1]=arr[i];
				}
				arr[i+1]=a;
			}
			if (jstack == 0) break;
			ir=istack[jstack--];			/* Pop stack and begin a new round of parti- */
			l=istack[jstack--];				/*	tioning. */
		} else {
			k=(l+ir) >> 1;
			SWAP(arr[k],arr[l+1])
			if (arr[l+1] > arr[ir]) {
				SWAP(arr[l+1],arr[ir])
			}
			if (arr[l] > arr[ir]) {
				SWAP(arr[l],arr[ir])
			}
			if (arr[l+1] > arr[l]) {
				SWAP(arr[l+1],arr[l])
			}
			i=l+1;							/* Initialize pointers for partitioning. */
			j=ir;
			a=arr[l];						/* Partitioning element. */
			for (;;) {						/* Beginning of innermost loop. */
				do i++; while (arr[i] < a);	/* Scan up to find element > a. */
				do j--; while (arr[j] > a);	/* Scan down to find element < a. */
				if (j < i) break;			/* Pointers crossed. Partitioning complete. */
				SWAP(arr[i],arr[j]);		/* Exchange elements. */
			}								/* End of innermost loop. */
			arr[l]=arr[j];					/* Insert partitioning element. */
			arr[j]=a;
			jstack += 2;
			/* Push pointers to larger subarray on stack, process smaller subarray immediately. */
			if (jstack > NSTACK) nrerror("NSTACK too small in sort.");
			if (ir-i+1 >= j-l) {
				istack[jstack]=ir;
				istack[jstack-1]=i;
				ir=j-1;
			} else {
				istack[jstack]=j-1;
				istack[jstack-1]=l;
				l=i;
			}
		}
	}
	free_ivector(istack,1,NSTACK);
}

#undef M
#undef NSTACK
#undef SWAP
#undef NRANSI






/**********************************************************************************
 Indexes an array arrin[1..n], i.e., outputs the array
 indx[1..n] such that arrin[indx[j]] is in ascending order
 for j = 1;2; : : :;N. The input quantities n and arrin
 are not changed.
**********************************************************************************/
void indexxd(int n, double *arrin, int *indx) {
	int l, j, ir, indxt, i;
	double q;

	/* initialize index vector with [1, 2, ..., n] */
	for (j = 1; j <= n; j++)
		indx[j] = j;

	l = (n >> 1) + 1;
	ir = n;

	for (;;) {
		if (l > 1) q = arrin[(indxt = indx[--l])];
		else {
			q = arrin[(indxt = indx[ir])];
			indx[ir] = indx[1];
			if (--ir == 1) {
				indx[1] = indxt;
				return;
			}
		}
		i = l;
		j = l << 1;
		while (j <= ir) {
			if (j < ir && arrin[indx[j]] < arrin[indx[j + 1]]) j++;
			if (q < arrin[indx[j]]) {
				indx[i] = indx[j];
				j += (i = j);
			}
			else j = ir + 1;
		}
		indx[i] = indxt;
	}
}






/***********************************************************************************
  Returns a normally distributed deviate with zero mean and unit variance
***********************************************************************************/
float gasdev(long *idum)
/* Returns a normally distributed deviate with zero mean and unit variance,
   using ran1(idum) as the source of uniform deviates.	*/

{
	float ran1(long *idum);
	static int iset=0;
	static float gset;
	float fac,rsq,v1,v2;

/* 	if (*idum < 0) iset=0;	*/
	if  (iset == 0) {					/* We don't have an extra deviate handy, so */
		do {
			v1=2.0*ran1(idum)-1.0;		/* pick two uniform numbers in the square ex- */
			v2=2.0*ran1(idum)-1.0;		/* 	tending from -1 to +1 in each direction, */
			rsq=v1*v1+v2*v2;			/* see if they are in the unit circle, */
		} while (rsq >= 1.0 || rsq == 0.0);	/* 	and if they are not, try again. */
		fac=sqrt(-2.0*log(rsq)/rsq);
		/* Now make the Box-Muller transformation to get two normal deviates. Return one and */
		/* save the other for next time. */
		gset=v1*fac;
		iset=1;					/* Set flag. */
		return v2*fac;
	} else {					/* We have an extra deviate handy, */
		iset=0;					/* so unset the flag, */
		return gset;			/* and return it. */
	}
}




// primary sorting:   A
// secondary sorting: B

#include <stdio.h>
#include "nrutil.h"
#include "mathfunctions.h"

int sort2(long n, double *A, double *B) {

    double *x, *y, similar;
    long i, j, index;
    long patches_x=1, patches_y=1;  // number of patches in x- and y-direction

/******************************************************************************/

    // create working vectors [1,n]
    x = dvector(1,n);
    y = dvector(1,n);

    for(i=1;i<=n;i++) {
        x[i] = A[i-1];
        y[i] = B[i-1];
    }


/******************************************************************************/

    // sort both vectors according to x
    for(i=1;i<=n;i++) {
        dmin(n-i+1,&x[i],&index);     // first element index = 0
        dshift(&x[i], &x[i+index]);
        dshift(&y[i], &y[i+index]);
    }


    similar   = x[1];
    patches_y = 1;
    for(i=2;i<=n;i++) {
        if (x[i] == similar) patches_y++;
        else break;
    }

    patches_x = n/patches_y;

    printf("sort2: found %ldx%ld matrix (collumns x rows)\n", patches_x, patches_y);


    // sort y-vector
    for (i=1;i<=patches_x;i++) {
        for (j=1;j<=patches_y;j++) {
            dmin(patches_y-j+1,&y[(i-1)*patches_y+j],&index);
            dshift(&y[(i-1)*patches_y+j], &y[(i-1)*patches_y+j+index]);
        }
    }


    for(i=1;i<=n;i++) {
        A[i-1] = x[i];
        B[i-1] = y[i];
    }



    /*printf("sorted vectors x and y\n");
    for(i=1;i<=n;i++) {
        printf("%lf  %lf\n", x[i], y[i]);
    }*/

    free_dvector(x,1,n);
    free_dvector(y,1,n);

    return(0);

}

